suppressPackageStartupMessages({
library(SummarizedExperiment)
library(SEtools)
library(edgeR)
library(DT)
library(pheatmap)
library(plotly)
library(dplyr)
library(sva)
library(ggrepel)
})
source("Functions/EDC_Functions.R")
source("Functions/CriFormatted.R")
load("Data/AllSEcorrected.RData",verbose=T)## Loading objects:
## SEs
## DEAs
load("Data/HormonalGenes/HormonalGenes.RData",verbose = T)## Loading objects:
## Thyroid
## Androgen
## Estrogen
## Corticoid
## Progesterone
## PPAR
## Retinoic
load("Data/DEAfetSingleCompunds.RData", verbose=TRUE)## Loading objects:
## res.fetT3
## res.fetBPA
load("Data/DEAorgSingleCompunds.RData", verbose=TRUE)## Loading objects:
## res.orgT3
## res.orgBPA
fetal <- SEs$chronic.fetal[,which(SEs$chronic.fetal$EXPO=="CNT" | SEs$chronic.fetal$EXPO=="1X" | SEs$chronic.fetal$EXPO=="1000X" | SEs$chronic.fetal$EXPO=="T3"| SEs$chronic.fetal$EXPO=="BPA0.04X")]
degs.fe <- row.names(DEAs$chronic.fetal)[which((abs(DEAs$chronic.fetal$logFC.EXPO1X)>0.5 | abs(DEAs$chronic.fetal$logFC.EXPO1000X)>0.5) & DEAs$chronic.fetal$FDR<=0.05 & DEAs$chronic.fetal$logCPM>0)]
sehm(fetal, assayName = "corrected", degs.fe, do.scale = T, show_rownames = F, anno_columns = c("EXPO2","MixtureBatch"), main="MixNDEGs")sehm(fetal[,which(SEs$chronic.fetal$EXPO=="CNT" | SEs$chronic.fetal$EXPO=="T3")], assayName = "corrected", Thyroid, do.scale = T, show_rownames = F, anno_columns = c("EXPO2","MixtureBatch"), main="T3 target genes")degs.fetT3 <- row.names(res.fetT3)[which(abs(res.fetT3$logFC)> 0.5 & res.fetT3$FDR<=0.05 & res.fetT3$logCPM > 0)]
sehm(fetal, assayName = "corrected", degs.fetT3, do.scale = T, show_rownames = F, anno_columns = c("EXPO2","MixtureBatch"), main="T3DEGs")T3 <- res.fetT3
T3$genes <- row.names(T3)
MixN <- DEAs$chronic.fetal
MixN$genes <- row.names(MixN)
MixN$logFC <- (MixN$logFC.EXPO1X + MixN$logFC.EXPO1000X)/2
Comp <- compareResultsFCNew(T3, MixN, FDRth=0.05, FCth=2^0.5, FDRceil=1e-10, logCPMth=0, FCceil=5,
title='Fetal: Comparing T3 to MixN (mean FC)', geneLabel=TRUE, topLab=-30)
Comp$Scatter#ggsave(Comp$Scatter, filename='Data/SingleCompounds/ScatterFetal_MixN_T3.pdf', width=10, height=10)
#ggsave(Comp$Scatter, filename='Data/SingleCompounds/ScatterFetal_MixN_T3.png', width=10, height=10)# Check overlap between MixN and Thyroid
table(Comp$AllRes$Status)##
## NotSignificant SignificantA SignificantA&B SignificantB
## 12488 374 17 44
## Fisher Test Significant
SigBoth <- dplyr::filter(Comp$Significant, Comp$Significant$Status=='SignificantA&B')
Contingency <- matrix(c(dim(dplyr::filter(SigBoth, log2FCA < 0 & log2FCB < 0 ))[1],
dim(dplyr::filter(SigBoth, log2FCA >= 0 & log2FCB < 0 ))[1],
dim(dplyr::filter(SigBoth, log2FCA < 0 & log2FCB >= 0 ))[1],
dim(dplyr::filter(SigBoth, log2FCA >= 0 & log2FCB >= 0 ))[1]),
nrow=2, dimnames=list(c('DownA', 'UpA'), c('DownB', 'UpB')))
Contingency## DownB UpB
## DownA 11 3
## UpA 1 2
fisher.test(Contingency)##
## Fisher's Exact Test for Count Data
##
## data: Contingency
## p-value = 0.1912
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.2518663 468.7960007
## sample estimates:
## odds ratio
## 6.321431
degs.fetBPA <- row.names(res.fetBPA)[which(abs(res.fetBPA$logFC)> 0.5 & res.fetBPA$FDR<=0.05 & res.fetBPA$logCPM > 0)]
sehm(fetal, assayName = "corrected", degs.fetBPA, do.scale = T, show_rownames = F, anno_columns = c("EXPO2","MixtureBatch"), main="BPADEGs")BPA <- res.fetBPA
BPA$genes <- row.names(BPA)
MixN <- DEAs$chronic.fetal
MixN$genes <- row.names(MixN)
MixN$logFC <- MixN$logFC.EXPO1X
Comp <- compareResultsFCNew(BPA, MixN, FDRth=0.05, FCth=2^0.5, FDRceil=1e-10, logCPMth=0, FCceil=4,
title='Fetal: Comparing BPA to MixN 1X', geneLabel=TRUE, topLab=-30)
Comp$Scatter#ggsave(Comp$Scatter, filename='Data/SingleCompounds/ScatterFetal_MixN_BPA.pdf', width=10, height=10)
#ggsave(Comp$Scatter, filename='Data/SingleCompounds/ScatterFetal_MixN_BPA.png', width=10, height=10)Organoids <- SEs$chronic.org[,which(SEs$chronic.org$EXPO=="CNT"|SEs$chronic.org$EXPO=="1X" | SEs$chronic.org$EXPO=="1000X" | SEs$chronic.org$EXPO=="T3"| SEs$chronic.org$EXPO=="BPA0.04X")]
degs.org <- row.names(DEAs$chronic.org)[which((abs(DEAs$chronic.org$logFC.EXPO1X)>0.5 | abs(DEAs$chronic.org$logFC.EXPO1000X)>0.5) & DEAs$chronic.org$FDR<=0.05 & DEAs$chronic.org$logCPM>0)]
sehm(Organoids, assayName = "corrected", degs.org, do.scale = T, show_rownames = F, anno_columns = c("EXPO2","MixtureBatch"), main="MixNDEGs")sehm(Organoids[,which(SEs$chronic.org$EXPO=="CNT" | SEs$chronic.org$EXPO=="T3")], assayName = "corrected", Thyroid, do.scale = T, show_rownames = F, anno_columns = c("EXPO2","MixtureBatch"), main="T3 target genes")degs.orgT3 <- row.names(res.orgT3)[which(abs(res.orgT3$logFC)> 0.5 & res.orgT3$FDR<=0.05 & res.orgT3$logCPM > 0)]
sehm(Organoids, assayName = "corrected", degs.orgT3, do.scale = T, show_rownames = F, anno_columns = c("EXPO2","MixtureBatch"), main="T3DEGs")T3 <- res.orgT3
T3$genes <- row.names(T3)
MixN <- DEAs$chronic.org
MixN$genes <- row.names(MixN)
MixN$logFC <- (MixN$logFC.EXPO1X + MixN$logFC.EXPO1000X)/2
Comp <- compareResultsFCNew(T3, MixN, FDRth=0.05, FCth=2^0.5, FDRceil=1e-10, logCPMth=0, FCceil=6,
title='Comparing T3 to MixN (mean FC)', geneLabel=TRUE, topLab=-40)
Comp$Scatter#ggsave(Comp$Scatter, filename='Data/SingleCompounds/ScatterOrg_MixN_T3.pdf', width=10, height=10)
#ggsave(Comp$Scatter, filename='Data/SingleCompounds/ScatterOrg_MixN_T3.png', width=10, height=10)# Check overlap between MixN and Thyroid
table(Comp$AllRes$Status)##
## NotSignificant SignificantA SignificantA&B SignificantB
## 12331 2212 60 286
## Contingency table on shared genes
SigBoth <- dplyr::filter(Comp$Significant, Comp$Significant$Status=='SignificantA&B')
ContingencyII <- matrix(c(dim(dplyr::filter(SigBoth, log2FCA < 0 & log2FCB < 0 ))[1],
dim(dplyr::filter(SigBoth, log2FCA >= 0 & log2FCB < 0 ))[1],
dim(dplyr::filter(SigBoth, log2FCA < 0 & log2FCB >= 0 ))[1],
dim(dplyr::filter(SigBoth, log2FCA >= 0 & log2FCB >= 0 ))[1]),
nrow=2, dimnames=list(c('DownA', 'UpA'), c('DownB', 'UpB')))
ContingencyII## DownB UpB
## DownA 14 36
## UpA 10 0
prop.table(ContingencyII)## DownB UpB
## DownA 0.2333333 0.6
## UpA 0.1666667 0.0
fisher.test(ContingencyII)##
## Fisher's Exact Test for Count Data
##
## data: ContingencyII
## p-value = 2.601e-05
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.0000000 0.2092738
## sample estimates:
## odds ratio
## 0
dplyr::filter(SigBoth, log2FCA < 0 & log2FCB > 0 ) %>% dplyr::pull(genes)## [1] "SP5" "COL9A2" "KRT8" "DDIT4L" "IGSF21" "NEUROG1"
## [7] "ID1" "HIST2H3A" "HIST2H3C" "CDC20B" "METRN" "ESM1"
## [13] "SPATA18" "RPE65" "SLC35F2" "TRIL" "HIST2H3D" "TRAPPC5"
## [19] "TMEM121" "S1PR3" "SALL3" "FEZF2" "PTHLH" "NDUFS7"
## [25] "TPH1" "SOD3" "ZNF503" "TGFB2" "PKMYT1" "DMRT3"
## [31] "COL9A3" "C4orf48" "CCDC85B" "APOE" "COL18A1" "C19orf24"
degs.orgBPA <- row.names(res.orgBPA)[which(abs(res.orgBPA$logFC)> 0.5 & res.orgBPA$FDR<=0.05 & res.orgBPA$logCPM > 0)]
sehm(Organoids, assayName = "corrected", degs.orgBPA, do.scale = T, show_rownames = F, anno_columns = c("EXPO2","MixtureBatch"), main="BPADEGs")BPA <- res.orgBPA
BPA$genes <- row.names(BPA)
MixN <- DEAs$chronic.org
MixN$genes <- row.names(MixN)
MixN$logFC <- MixN$logFC.EXPO1X
#Comp <- compareResultsFC(BPA, MixN, FDRth=0.05, FCth=2^0.5, FDRceil=1e-10, FCceil=8,logCPMth = 0, title='Comparing BPA to MixN 1X', geneLabel=TRUE)
#Comp$Scatter
Comp <- compareResultsFCNew(BPA, MixN, FDRth=0.05, FCth=2^0.5, FDRceil=1e-10, logCPMth=0, FCceil=4,
title='Organoids: Comparing BPA to MixN 1X', geneLabel=TRUE, topLab=-30)
Comp$Scatter#ggsave(Comp$Scatter, filename='Data/SingleCompounds/ScatterOrg_MixN_BPA.pdf', width=10, height=10)
#ggsave(Comp$Scatter, filename='Data/SingleCompounds/ScatterOrg_MixN_BPA.png', width=10, height=10)# Check overlap between MixN and BPA
table(Comp$AllRes$Status)##
## NotSignificant SignificantA SignificantA&B SignificantB
## 13770 821 50 248
## Contingency table on shared genes
SigBoth <- dplyr::filter(Comp$Significant, Comp$Significant$Status=='SignificantA&B')
ContingencyII <- matrix(c(dim(dplyr::filter(SigBoth, log2FCA < 0 & log2FCB < 0 ))[1],
dim(dplyr::filter(SigBoth, log2FCA >= 0 & log2FCB < 0 ))[1],
dim(dplyr::filter(SigBoth, log2FCA < 0 & log2FCB >= 0 ))[1],
dim(dplyr::filter(SigBoth, log2FCA >= 0 & log2FCB >= 0 ))[1]),
nrow=2, dimnames=list(c('DownA', 'UpA'), c('DownB', 'UpB')))
ContingencyII## DownB UpB
## DownA 32 0
## UpA 5 13
prop.table(ContingencyII)## DownB UpB
## DownA 0.64 0.00
## UpA 0.10 0.26
fisher.test(ContingencyII)##
## Fisher's Exact Test for Count Data
##
## data: ContingencyII
## p-value = 2.414e-08
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 12.94172 Inf
## sample estimates:
## odds ratio
## Inf
N.B. Overlap analyses with the function below are not equivalent to the scatterplots, because here the selection of DEGs for MixN is done with an ‘or’ approach among the two concentrations for the FC (while above we use MixN 1x FC for BPA and the mean between the two concentrations for T3).
library(overlapper)##
## Attaching package: 'overlapper'
## The following object is masked _by_ '.GlobalEnv':
##
## overlap.prob
m1 <- overlapper::multintersect(ll=list(fetalDEGs=degs.fe, fetalBPA=degs.fetBPA, fetalT3=degs.fetT3), universe=row.names(DEAs$chronic.fetal), two.tailed = F)
dotplot.multintersect(m1, sizeRange = c(0,15), th=0.05)## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_text).
m2 <- overlapper::multintersect(ll=list(orgDEGs=degs.org, organoidBPA=degs.orgBPA, organoidT3=degs.orgT3), universe=row.names(DEAs$chronic.org),two.tailed = F)
dotplot.multintersect(m2, sizeRange = c(0,15), th=0.05)## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_text).